%% Clean up
clear variables

%% General
type    = 'SA';
aux.n_m = 10^3;
aux.n   = 10^3;

%% Paths
filename    = 'chains.xlsx';
datapath    = 'Output\';
path        = 'parameters_SA\';

%% Exo parameters and moments
param.shock     = 0.01;% Size of shock
param.r         = log(1.05)/12;% interest rate
param.alpha     = 0.64;% 
param.epsilon   = 0.5;% matching efficiency

%% Moments
param.u_mom             = 0.06;% unemployment rate moment
param.lambda_mom        = 0.25;% job offer arrival rate
param.lambda            = param.lambda_mom;
param.zeta              = 0;% exo separation rate
param.job_to_job_mom    = 0.032;% job to job rate
param.mom_deudAPL       = -3.6;% cyclicality of EU rate
param_ex                = param;
disc_lis                = nan(2,1);

for i = 1:2
    param = param_ex;
    if i == 1
        name = type;
        theta0=[ 0.28, 0.5, 40];   
        param.inac_mom          = 0.225;% 12m employment weighted inaction rate
    else
        name = ['OJS_',type];
        theta0=[ 0.15, 0.1];
        param.inac_mom          = 0;% 12m employment weighted inaction rate
    end

    theta0 = fminsearch(@ (x) moments_SA( x, param, aux, 1), theta0, optimset('display', 'iter', 'MaxIter',50));
    theta0 = fminsearch(@ (x) moments_SA( x, param, aux, 1), theta0, optimset('display', 'iter'));
    disc_lis(i) = moments_SA( theta0, param, aux, 1);

    param = moments_SA( theta0, param,aux, 0);
    param_ini = param;

    x0 = 0.991;
    x0 = fzero(@ (x) root_w_SA( x, 1, param, param_ini, aux), [x0^2,1], optimset('display','iter'));
    x0 = fzero(@ (x) root_w_SA( x, 1, param, param_ini, aux), [x0-0.001, x0+0.001], optimset('display','iter'));
    param_end = root_w_SA( x0,0,param,param_ini, aux);
    param_ini.epsilon = param_end.epsilon;

    save( [path,'param_ini_',name,'.mat'], 'param_ini');
    save( [path,'param_end_',name,'.mat'], 'param_end');


    dlnu_dlnAPL   = log(u_fun(param_ini)./u_fun(param_end))./log(mean_fun( aux.n, 'm', param_ini)./mean_fun( aux.n, 'm', param_end));
    dlnlam_dlnAPL = log(param_ini.lambda./param_end.lambda)./log(mean_fun( aux.n, 'm', param_ini)./mean_fun( aux.n, 'm', param_end));
    dlnEU_dlnAPL  = log(u_fun(param_ini).*param_ini.lambda/(1-u_fun(param_ini)) /...
        (u_fun(param_end).*param_end.lambda/(1-u_fun(param_end))))./log(mean_fun( aux.n, 'm', param_ini)./mean_fun( aux.n, 'm', param_end));

    m_ini = linspace( param_ini.m_l, param_ini.m_u, aux.n)';
    m_end = linspace( param_end.m_l, param_end.m_u, aux.n)';

    dlnEE_dlnAPL= log( jtj_num( delta_fun( m_ini, param_ini), q_fun( m_ini, param_ini))/ ...
        jtj_num( delta_fun( m_end, param_end), q_fun( m_end, param_end)))./log(mean_fun( aux.n, 'm', param_ini)./mean_fun( aux.n, 'm', param_end));

    param_ini.A = param.lambda_mom;% assign matching efficiency
    param_end.A = param_ini.A;

    dlnV_dlnAPL=log(v_fun(u_fun(param_ini),param_ini.lambda,param_ini)./v_fun(u_fun(param_end),param_end.lambda,param_end))./log(mean_fun(aux.n,'m',param_ini)./mean_fun(aux.n,'m',param_end));

    T = table( dlnu_dlnAPL, dlnV_dlnAPL, dlnlam_dlnAPL, dlnEU_dlnAPL, dlnEE_dlnAPL);
    
    if param.K == 0
        writetable(T,[datapath,filename],'Sheet','Table C','Range','A4')
    else 
        writetable(T,[datapath,filename],'Sheet','Table C','Range','A1')    
    end

    clear param_ini param param_end
end